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VIII. 


I. INTRODUCTION 


Numerous studies in recent years have demonstrated the usefulness of 
the PNS equations In the calculation of a definable class of supersonic 
flows. If the inviscid region of the flow Is supersonic and there is no 
streamwise separation of the flow, the equations of motion can be accurate- 
ly modelled by a mixed set of hyperbolic-parabolic equations (the PNS 
equations). These equations can be solved much more efficiently than the 
complete Navier-Stokes equations since the solution can be marched in 
space rather than title. 

Various versions of PNS equations have been successfully employed. 

One of the earliest studies involving the PNS equations was performed by 
Rudman and Rubin (1) in 1968. Rudman and Rubin applied a series expansion 
technique to the steady Navier-Stokes equations and by eliminating higher- 
order terms produced a system of strictly parabolic Navier-Stokes equa- 
tions. A less formal approach was taken by Lubard and Helliwell (2) in 
which streamwise viscous stresses were assumed small in comparison with 
the normal viscous stresses. Thus, the Lubard-Helliwell PNS system is 
derived' by dropping viscous terms containing partial derivatives in the 
streamwise direction. The retention of the pressure gradient in the 
streamwise momentum equation of this system is the most significant dif- 
ference between the Lubard-Helliwell PNS equations and the Rudman-Rubin 
PNS equations. The absence of this pressure gradient term allows for 
stable space marching but may lead to inaccuracies in flowflelds containing 
moderate streamwise pressure gradients. In this investigation, the more 



2 


common Lubard-Helliwell formulation is employed with the pressure gradient 
term being treated in a manner described in the next section. 

The PNS equations have been integrated using a variety of finite- 
difference schemes. Because of its ease of implementation, a simple 
explicit scheme was used by Rudman and Rubin (1) for their calculations 
of the merged layer region near sharp leading edges in hypersonic vis- 
cous flow. In an effort to obtain solutions farther downstream, Rubin and 
Lin (3) proposed a predictor-corrector, semi-implicit, multiple-iteration 
scheme. Due to the larger allowable marching step size, this scheme was 
found to require an order-of magnitude less computer time to perform the 
same calculations than the explicit scheme. In their investigation of 
hypersonic viscous flow over cones at high angle of attack, Lubard and 
Helliwell (2) used an implicit differencing of the equations with a 
Newton-Raphson iteration technique to solve the resulting systems of non- 
linear algebraic equations. In the late 1970s, noniterative, implicit, 
approximate-factorization schemes were developed by Vigneron et al. (4) 
and Schlff and Steger (5). These schemes were based on a class of ADI 
schemes developed by Lindemuth and Killeen (6), McDonald and Briley (7), 
and Beam and Warming (8). Though they require the inversion of block 
tridiagonal systems of linear algebraic equations in the calculation of 
flow properties at each step, these schemes were found to be more compu- 
tationally efficient than the iterative schemes previously used, and they 
are the schemes most commonly employed in PNS calculations today. 

In 1981, MacCormack (9) proposed an implicit scheme which requires 
only the inversion of block bidiagonal systems rather than block tri- 
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diagonal systems, thus yielding a savings in computer time and storage 
requirements. This scheme was designed to solve time dependent equations 
such as the complete Navier-Stokes equations. It is based on MacCormack,' s 
well-proven second-order accurate explicit predictor-corrector method (10) 
but adds an implicit procedure in the predictor-corrector sequence for 
points in the flow at which the local CFL number exceeds the stability 
limit. The method has been applied to two-dimensional internal supersonic 
flows (11, 12), two-dimensional external flows (13), external axisymmetric 
flows (1 h), quasi-one-dlmensional flows (15), three-dimensional flow over 
a biconic with compression flap (16), as well as three-dimensional blunt 
body flows. In each of these cases, the scheme was applied to either the 
complete, or thin layer forms of the unsteady, Navier-Stokes equations 
as well as the viscous shock layer equations. 

In the present work, the implicit MacCormack scheme has been modified 
to solve the parabolized Navie7.‘-Stokes equations. This report describes 
the resulting finite-difference algorithm and presents computational 
results for two laminar test cases. Results for the case of a flat plate 
boundary layer are compared with those obtained using the conventional 
Beam-Warming scheme as well as those obtained from a boundary -layer code. 
In a more severe test of the method, the hypersonic flow past a 15° com- 
pression corner has been computed. For this case, a global iteration on 
the pressure field of the form developed by Raklch (17) was applied in 
conjunction with both the implicit MacCormack scheme and the Beam-Warming 
scheme. Using an iteration of this type, it is possible to include 
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Influences from downstream which are otherwise neglected in a parab- 
olized Navier-Stokes calculation. The computed results are compared 
with available experimental data and a numerical solution of the com- 
plete. Navier-Stokes equations. » 


* 
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II. GOVERNING EQUATIONS 


A. Navler-Stokes Equations 

The equations describing the planar flow of a Newtonian fluid are 
the two-dimensional, unsteady Navier-Stokes equations. These ca"» >-«? 
written in nondimenslonal strong-conservation-law form in Cartesian coor'-* 


dinates as 
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The equations have been nondlmenslonallzed (dimensional quantities are 
denoted by a tilde) in the following manner: 
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P = P/P« V ~ 

where L is the reference length of unity. 

The Reynolds number (Rs^) appearing In the viscous terms is given by 


Re 


00 



The coefficient of thermal conductivity has been replaced by assuming a 

constant Prandtl number and the coefficient of viscosity is calculated 

using Sutherland's equation 

/9 / l + 110.4 k/t\ 
y " T \T + 110.4 K/T m ) 

Finally, the system is closed using the perfect gas equation of state 
which in nondimens ional form becomes 


p = pT/yM* 
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B. Coordinate Transformation 
A transformation of the spatial coordinates of the form 
5 » x 


n - n(x,y) 


Is applied to Eq. (1) so that the equation may be differenced on a uniform 
computational mesh. The resulting equation In strong-conservation-law 
form is 


A 

bjj 

3 t 


+ 
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( 2 ) 


where 

0 = u/j 

E = 7 (E - E.) 

J V 

F = i [n x (E - E v ) + n y (F - F v )] 

and J Is the Jacobian of the transformation given by 

t - 3(€.n) 

J " 3<x,y) 


Derivatives in the viscous vectors are transformed using the chain 


rule: 


3 F 3 _ _ 3 _ 

3x “ S 3^ “ H 

3_ _ , n 3 

3y " n x 3x n y 3y 


C. Parabolizing Assumptions 

The equations are "parabolized" to permit stable marching in space 
by making the following assumptions : 
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1. Steady flow. 

2. Strenmwlse viscous derivatives are negligible in comparison with 
normal viscous derivatives. This approximation is valid for flows with 
high Reynolds numbers. 

The following system of PNS equations Is obtained as a result of 
these assumptions 


3E x ' aF n 
+ ° 


(3) 


where 


E - ~ E 

F - 7 In (E - EJ + n,.(F - F ) 3 
J x v y v 

and E v and F v now contain no £ derivative terms. The PNS equations are 
a mixed set of hyperbolic-parabolic equations In the streamwlse direction 
5 provided that the invlscld flow Is supersonic, the streamwise velocity 
component Is everywhere greater than zero, and the streamwlse pressure 
gradient term in the streamwise momentum equation is either omitted or 
the "departure behavior" Is suppressed by a suitable technique. 


D. Streamwise Pressure Gradient 

The presence of the streamwlse pressure gradient term in the stream- 
wise momentum equation permits information to be propagated upstream 
through subsonic portions of the flowfield such as a boundary layer. As 
a consequence, a space-marching method of solution Is not well-posed 
and in some cases exponentially growing solutions (departure solutions) 
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arc encountered, A number of different techniques have been proposed 
to eliminate this difficulty. For this study, the method proposed by 
Vigneron et al. (4) is used. 

The "Vigneron technique" involves splitting the E vector into two 
parts, 

E * E* + P 

where 

E* « “ [pu, pu 2 f uip , puv , (e t + p)u] T 
? - j [0, (1 - w)p, 0, 0] T 


The E* vector now replaces E In the numerical scheme and Is treated 
as a source term which is evaluated in the supersonic region. The final 
form of the governing equations becomes 


BE* , BF 
H 3n 
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An eigenvalue analysis shows that this system will be hyperbolic-parabolic 
if 

2 
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where M- is the local streamwise Mach number. Since this relation was 
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determined using a linear analysis, a safety factor 0 is applied ant? 03 
is calculated from 
OYM 2 

03 = § 2 
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or Is set equal to one when the computed U) is greater than one. 


XIX, NUMERICAL SOLUTION OF EQUATIONS 


A. Numerical Scheme 

The numerical scheme used In the present study to Integrate the FNS 
equations Is an adaptation of the method proposed by MacCormack (9) In 
1981. MacCormack demonstrated the method's usefulness for solving the 
full unsteady Navier-Stokes equations In an application to a two-dimen- 
sional shock-boundary-aayer interaction problem. The method Is implicit 
In nature and thus allows a much larger marching step size than explicit 
methods. In addition, the method possesses three advantages over con- 
ventional fully implicit methods. First, the method uses two-point, one- 
sided differences In the implicit part of the algorithm. Thus, block 
bldiagonal systems of algebraic equations result which are significantly 
less i ‘4 <H,ly to invert than the block tridiagonal systems found In conven- 
tional methods. Second, the method employs the lnviscid Jacoblans and 
corrects them using representative viscous terms added to the Eulerlan 
eigenvalues. This maintains stability while avoiding the expensive cal- 
culation of the complete viscous Jacoblans. Finally, the algorithm 
allows the Implicit step to be skipped in regions where the explicit 
stability restriction is satisfied, as in the region away from the bound- 
ary layer where mesh spacing is large. The method is stable for unbounded 

y At 

At and second-order accurate under the condition that ~ 9 n 

P min (Ax , Ay ; 

remains bounded as a 2-D Cartesian mesh Is refined. 

In adapting this scheme for use in solving Eq. (4) the procedure 
outlined in Ref. 9 is followed. Equation (4) is differentiated with 


respect, to E, to obtain the equation 
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which governs the propagation of the local changes, In Eq. (5), 

C is the Jacobian 3 F/3E*. If the ^-derivative i s Implicitly differenced, 
the following equation results 


/. , ..BC'l/lw"* 1 /3E \ n / 3 P \ n+1 ... 

( I + i5 3— )\w) ■ m) -\w <6) 

where the dot in the equation indicates that the derivative also operates 
on the factor to the right. Next, define 



and 



so that Eq. (6) becomes 

( I + 45 |£l ) SE* n+1 = AE n -AP n+1 


(7) 


At this point, it is observed that proceeding to difference this 
equation in the same manner that MacCormack used for the unsteady equa- 
tions will yield a scheme which requires the modal analysis (i.e., finding 
the eigenvalues and eigenvectors) of the C matrix. In addition, the 
difficult task of decoding the E* vector will be encountered at each step 
in order to obtain the flow variables. 

This problem consists of solving the nonlinear set of algebraic 


equations 
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for the dependent variables, p, u, v, and p. This system can be reduced 
to a single quadratic equation for u-velocity, with one root of this equa- 
tion being associated with a supersonic flow and the other one a subsonic 
flow. Solution is straightforward when it is known a priori that the flow 
is supersonic everywhere as when integrating the Euler equations. However, 
for viscous problems, the solution is more difficult due to the presence 
of the subsonic region of the boundary layer. For a boundary -layer calcu- 
lation, the E* vector may be readily decoded at every point normal to the 
wall except at the point closest to the sonic line. It is possible to 
determine this point in most cases by linearly interpolating the Mach num- 
ber between the grid points on either side of the unknown point. The E* 
vector for this point is then decoded by choosing the branch of the solu- 
tion yielding a Mach number closest to the interpolated line. However, 
this is an undesirable method from the standpoint of computer time and 
code robustness. 

In order to avoid this procedure, as well as the modal analysis of 
the C matrix, Eq. (7) is modified using the linearization 

< S I *," +1 = A 5 (| li ) n (||) n+1 + 0 ( A e , 2 
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and substituting Eq. (9) into Eq. (7), we have 

(a + AZ • 'j 6fj n+1 = Al n - AP n+1 (10) 

where 

B = i| 

3U 

The equation now contains the simpler Jacobians, A and B, and the easily 
decoded U vector. Additional justification for this change in the depen- 
dent variable is provided by consideration of the case of a more general 
transformation of the streamwise coordinate. For example, a transforma- 
tion of the form Z, - £(x, y) allows the marching planes to be of general 
shape and orientation, which is desirable for marching along bodies with 
large surface slope. However, this transformation also yields an E* 
vector which is virtually impossible to decode, thus making a change in 
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marching vector from E* to U essential. The disadvantages of the present 
formulation will be discussed in a subsequent paragraph. 

Differencing Eq. (10) in a manner consistent with MncCormack’s 
original explicit scheme (10) yields the following implicit predictor- 
corrector algorithm for the numerical integration of Eq. (4). 


Predictor: 


* An 


(lla) 


( A j ' A5 l B l , )' sC j +1 - aE ” 
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(lib) 
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Corrector: 


A F n+1 


(12a) 
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(12b) 
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The differencing operators, A + /An and A_/An, are defined by 


A. z z . . _ - z 
+ i+l 


A z z . - z . 
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The right-hand sides of Eqs. (11a) and (12a) are calculated as In the 
explicit scheme. That Is, one-sided differences are used for the convec- 
tive terms and central differences for the viscous terms. In the present 
code, the differencing permutation can be reversed at each marching step 
depending on the boundary condition used (see subsection C) . 

The matrix B in Eq. (10) has been replaced with the related matrix 
| B ( in the implicit steps of the algorithm. This substitution is required 
in order for the block bidiagonal systems to be inverted numerically in 
a stable manner. The matrix |b| is defined by 

l»l = S’ 1 D B S y (13) 


where is the matrix whose rows are the left eigenvectors of the in- 

viscid Jacobian, B, = 3F./9U (see the Appendix). D„ is the diagonal matrix 
1 i B 

whose elements d^ are defined by 
d B ± = K^J + VISCOR) 

where A^ is the ith eigenvalue of B^ (Appendix). The viscous correction, 
VISCOR, is related to the eigenvalues of the viscous Jacobian SF^/SU and 
is given by 


VISCOR = |^ 
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The coefficient K is determined by 
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where Re^ = Re^puAo/ (yhy) and a is a safety factor usually set equal to 


0.5. Also, is the largest step size which satisfies the CFL 

condition 

(A5 >cfx i r 1 

c 

where is the maximum eigenvalue of 3F^/3E*. Thus, when the step size 
is such that the explicit scheme is locally stable, K is set to zero and 
the Implicit step reduces to 

[A ]6U. «* AE, - AP. (14) 

y J 1 3 

This matrix equation reveals a disadvantage of marching with the U 
vector inasmuch as it is necessary to invert a 4 x 4 matrix at each point 
including those a which the explicit scheme is stable. The other draw- 
back of this formulation is that the main diagonal blocks of the coeffi- 
cient matrices which result cannot be easily diagonalized as MacCormack 
has done. A lower-upper decomposition is used here to invert the main 
diagonal blocks. It is believed, however, that the computer time saved 
by diagonalizing the coefficient matrices would be spent in decoding the 
E* vector. In order to remedy the first disadvantage, steps b) and c) 
of Eqs . (11) and (12) may be replaced with 

Predictor: 

E *?* 1 = E* n + AE? - AP? 

3 3 3 3 


Corrector: 


( 15 ) 


17 


in regions of the flow where flow properties and grid spacing are such 
that the marching step size satisfies the stability restriction of the 
explicit scheme. Generally, the explicit scheme will not be stable in 
the subsonic region of the boundary layer so that decoding is not a prob- 
lem unless a more general transformation has been applied. 

The implicit smoothing proposed by MacCormack was found to be un- 
necessary for the test cases of the present study. However, when cap- 
turing shocks some explicit damping is required. The pressure smoothing 
term of Ref. 18 was found to be sufficient to allow stable space marching. 
This scheme Involves adding a term proportional to the quantity 
(Ar))^ [p^lu so that the smoothing effect is most pronounced in regions 
of high pressure gradient changes. In the present code, the term is 
centrally differenced as 

aE * - c Ip hi - fr.1 + p i- il ( V - + Vl> 

p 1+1 + *p, + P,-! 


j 


where C is proportional to the marching step size. This term is fourth 

X 

order in Af| so that little effect on the viscous forces is observed. 


B. Computational Grid 

For viscous flow cases, it is desirable to attain as much resolution 
of the boundary layer as possible in order that the viscous stresses may 
be accurately modelled. This is accomplished in the present study by 
clustering the grid points normal to the wall using the following clus- 


tering function 
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z(n) - 


/ 3±1 \ 1 *“ n 

3+1 - (3-1) \ 3-1 / 



where z(r|) becomes equal to one when n » 1 and Is zero vhen r) « 0. The 
points become more tightly clustered as g approaches 1. 

The physical distance, y, is then obtained from 

y(?,n) - y 0 (5) + 6(0 *(n) 

where 6(0 is the distance from the wall to the freestream edge of the 
grid and y (O is the value of y at the wall. 

The transformation used is analytic and the metrics could therefore 
be determined analytically. However, experience has shown that it is 
desirable to calculate the terms n x and r] numerically in a manner con- 
sistent with the finite-difference scheme being used. For the HacCormack 
scheme, one-sided differences are used to evaluate the Jacobian of the 
transformation, J. The differencing follows that of the convective 
derivatives. The Jacobian can be calculated by the equation 



The metrics are then computed using the relations 



y 5 (n) - y 0? 


<0 + 6 ; <S) z(n) 


and 
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C. Boundary Conditions 

At the wall boundary, no-slip conditions were Imposed on the veloc- 
ities and the temperature was specified. In order to determine the wall 
pressure, a zero-gradient extrapolation was explicitly applied to the 
pressure at the end of each predictor and corrector sweep. 

At the freestream edge of the grid, zero-gradient extrapolations 
were explicitly applied to all of the dependent variables at the end of 
each predictor and corrector step. 

One of the weaknesses of the present scheme lies in the difficulty 
of obtaining a reliable Implicit boundary condition. That is, at the 
start of each predictor and corrector sweep a value for the vector, 

~~ | B | 6U, is required at the boundary. At boundaries where the flow 
properties and grid spacing are such that the explicit scheme is stable 
this presents no problem since |b| is zero. This Is generally the case 
at freestream boundaries. However, at solid wall boundaries In viscous 
flow, the grid spacing normal to the wall is generally small and implicit 
treatment is required. The approach generally taken is that described by 
MacCormack (9) in which the flux, (b|6u, crossing the boundary at the end 
of the predictor step is reinjected into the flow to start the corrector 
sweep. That is, 



where 
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[E] 


10 0 0 

0 10 0 

0 0-1 0 

0 0 0 0 


Of course, the use of this boundary condition requires that the differ- 
encing permutation for the right-hand side [Eqs, (11a) and (12a)] remain 
the same for each step such that the predictor step is always swept 
toward the wall and the corrector step away from the wall. 

The boundary condition described above was tested on both test cases 
and its performance was compared with that of the condition which simply 
sets | ^j* |B|6uj” + ^ equal to zero. The results will be discussed in the 
following section. 
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XV. NUMERICAL RESULTS 

In order to evaluate the present implicit method for solving the PNS 
equations, two test cases were computed, 

A. Test Case I 

For the first test case, the supersonic laminar flow over a flat 

plate (see Fig. 1) was computed. The freestream flow conditions for this 

case are 

M = 2.0 
00 

Re /L » 1.99 x 10 7 /m 
00 

^ = T “ 233 K 
00 w 

Pr = 0,72 

Two PNS codes were written to compute thi;' boundary-layer test case. 
The first code used the conventional Beam-Warming finite-difference scheme 
while the second employed the Implicit MacCormack scheme. The results 
from the two codes were compared with the results obtained from the com- 
pressible boundary-layer code of Pletcher (19). Thus, a fair evaluation 
of the suitability of using the Implicit MacCormack scheme to solve the 
PNS equations could be made. 

The initial conditions at x -• 1.54 were obtained from the boundary- 

layer code. To obtain the results presented here, the same equally 

-5 

spaced grid was used in all the calculations with Ay = 6.096 x 10 
Grid points were added to the top of the mesh as required by the growth 




of the boundary layer, This followed the procedure of the boundary-layer 
code and thus allowed a point by point comparison of the results. 


Profiles of tangential velocity and temperature at x * 4.57 are shown 
in Figs. 2 and 3, respectively. Figure 4 shows the streamwlse variation 
of skin friction coefficient as calculated from the formula 


C 


f 


M wall 9u 

Re w 5? 


Plotted In these three figures are results of the boundary-layer code, 
the results from the Beam-Warming PNS code, and the results from the 
Implicit MacCormack PNS code. The results are In excellent agreement. 
Though not shown, calculations were also performed using the explicit 
MacCormack scheme which proved to be unstable at a Courant number greater 
than 1.5, The PNS calculations were performed at a Courant number of 330. 

As indicated In the previous section for the Implicit MacCormack 
scheme, two different methods of treating the solid boundary were at- 

* 0 which 

caused the scheme to go unstable at Courant numbers larger than approxi- 
mately 409. The second method consists of reflection of the quantity 
|b|< 5U at the wall, This condition allows the new scheme to remain 

4 

stable at Courant numbers greater than 10 though calculated siiin fric- 
tion coefficients became inaccurate when the Courant number exceeded 

3 

5(10 ). The Beam-Warming code which employs straightforward implicit 

conditions at both boundaries was observed to be stable and to yield 

reasonable results for the skin friction coefficients at Courant numbers 
4 

of more than 10 . 


tempted. The first made use of the expression 
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Figure 2. Comparison of velocity profiles 
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Figure 3. Comparison of temperature profiles 
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Figure 4. Comparison of skin friction coefficients 




27 


The difference in accuracy of the two schemes was originally believed 
to be due to errors introduced by the approximation of the viscous terms 
on the left-hand side of the Implicit MacCormack algorithm. However, a 
comparison of the resultB obtained by using the two boundary condition 
procedures described above indicates the powerful effect of the boundary 
conditions on the stability and accuracy of the scheme. Observation of 
this effect leads one to conclude that the boundary condition procedure 
is likely to still be the dominant source of inaccuracies in implicit 
MacCormack scheme calculations on equally spaced grids. 

In addition to these calculations, some calculations were performed 

on stretched grids which varied in height with the streamwlse coordinate. 

Though the results are not included here, these calculations showed that, 

when using the implicit MacCormack scheme, the allowable Courant number 

is strongly dependent on the extent to which the grid is clustered in the 

3 

normal direction. In general, at Courant numbers larger than 10 , even 
slight amounts of grid stretching had catastrophic effects on the solu- 
tions. Calculations performed in the freestream with extrapolated 
boundary conditions (i.e., no wall boundary imposed) suggest that, with 
the new scheme, grid clustering has the effect of injecting nonphysical 
mass into the flow at the grid points. To eliminate this behavior, 

Hindman (20) emphasized the need to numerically solve the additional grid 
conservation-law equation, 
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Th^s equation must be solved with the same Integration scheme used for 
the transformed equations of motion when the equations are written in 
strong-conservation-law form. Satisfaction of this geometric conservation 
law is observed to correct the problem when using the explicit MacCormack 
scheme; however, extension to the implicit scheme is not straightforward 
and nonphysical behavior still persists in regions of the flow in which 
the implicit steps are required for numerical stability. These errors 
were also present in calculations using the Beam-Warming scheme on these 
grids; but because this scheme is centrally differenced, the errors are 
much less significant. 

Comparison was also made of the computer time required to perform 

-4 

the calculations. The Beam-Warming scheme required 9.629(10 ) sec/step/ 

grid point of CPU time on an NAS AS/6 computer. For the implicit Mac- 

_3 

Cormack scheme, 1.274(10 ) sec/step/grid point ware required. It must 

be pointed out, however, that the latter is a worst case value. That is, 
since the grid is evenly spaced, all points were calculated implicitly. 

A timing study was' performed to determine the reasons for the greater 
computer time of the implicit MacCormack scheme. As might be expected, 
a major contributor is the |b| matrix calculation which must be performed 
twice per marching step. Another significant factor is the calculation 
of the right-hand side terms which also must be computed twice per step. 
Thus, although the study showed that two block bidiagonal systems can be 
Inverted about 10% more quickly than one block tridiagonal system, some 
points must be calculate d explicitly for the implicit MacCormack scheme 
to be as efficient in solving the PNS equations as the Beam-Warming scheme. 
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A comparison of the computer storage requirements of the two schemes 
showed that the Beam-Warming scheme required 108 K bytes of storage while 
the implicit MacCormack scheme required only 76 K bytes. The reason for 
the lower requirements of the implicit MacCormack scheme is that block 
tridiagonal systems need not be stored for the entire grid. The block 
bldiagonal systems may be formed and inverted with one sweep while storing 
only two 4 x 4 matrices at a time. 

The two schemes were of about equal difficulty to code for the two- 
dimensional PNS equations. The difficulty of coding the viscous Jacobians 
in the Beam-Warming scheme is balanced by the coding of the |b| matrix 
in the MacCormack scheme. However, if a code exists which employs the 
explicit MacCormack scheme, the new scheme may be implemented by simply 
augmenting the existing algorithm with the Implicit steps. The Beam- 
Warming code of the presen?; study contained about 30 percent more FORTRAN 
source lines than the implicit MacCormack scheme mostly as a result of 
the block tridiagonal solver. 


B. Test Case II 

The second case computed was that of hypersonic laminar flow over a 

15° wedge. The flow conditions were chosen to correspond with one of the 

cases studied experimentally by Holden and Moselle (21) and numerically 

by Hung and MacCormack (18) using the complete Navier-Stokes equations. 

The flow conditions were 

M = 14.1 l = 0.439 m 

00 

T m = 72.2 K Pr = 0.72 

Re £ = 1.04 x 10 5 T w = 297 K 


The Reynolds number, Re^, Is the freestream Reynolds number based on the 
distance from the leading edge to the beginning of the ramp. This flow 
is supersonic in the inviscid region and exhibits no strearawise separation. 
Thus, stable space marching is allowed. The problem is illustrated sche- 
matically in Fig. 5. The grid used in the calculation is shown in Fig. 6 
with every other grid line omitted. Thirty grid points were spaced normal 
to the wall with a stretching parameter, $, of 1.08. The grid has an 
initial height of 0.139 t and the top of the grid is initially at an angle, 
<j> t , of 5° with respect to horizontal. At the beginning of the ramp, 
cj) to p changes to 15° to follow the rise of the wall. 

The Initial conditions at x = 0 were provided by specifying free- 
stream conditions everywhere except at the wall where no-slip conditions 

and constant wall temperature were imposed. The computation then pro- 

-3 

ceeded downstream with a step size of A B, = 3.05 x 10 and was terminated 
after 266 steps at x = 1.74 Z. About 13 seconds of CPU time on an NAS 
AS/6 computer were required for the calculation. This compares with the 
32 minutes of computer time on a CDC 7600 which were required by Hung and 
MacCormack to solve the complete Navier-Stokes equations. It should be 
noted that about 20 of the 30 points normal to the plate were computed 
explicitly by the present algorithm. 

Comparison of the wall pressures on the flat plate with those computed 
by Hung and MacCormack is shown in Fig. 7. Also presented in this figure 
are the theoretical results of the strong-interaction analysis of Bertram 
and Blackstock (22) . Good agreement is observed between the present 




Figure 6. Computational grid 
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results and those previously obtained by computational and theoretical 
methods. A small disagreement with the result of the experiment of Holden 
and Moselle is observed. The reason for the discrepancy is unknown. 

The pressure coefficients defined by C * p are compared with 

p w 

the previous computational and theoretical results in Fig. 8. Again, the 
present results compare well with those obtained by other methods. Be- 
cause of the single sweep marching, there is slight disagreement in the 
region near the beginning of the ramp since the flow upstream is not 
"warned" of the oncoming compression. Nevertheless, the results down- 
stream compare with the experiment as well as those obtained by the 
complete Navier-Stokes code. 

As noted by Hung and MacCormack, the intersection of the leading edge 
shock with the shock induced by the wedge results in a Type VI interfer- 
ence as classified by Edney (23). This interference produces an expansion 
fan which eventually Impinges on the wedge surface. However, this point 
of impingement falls downstream of the region computed in the present 
calculation, and therefore the peak of the computed pressure coefficient 
is expected to lie downstream of the final £ station. 

Heat transfer coefficients as calculated from 


_ ‘'w sec e 8T 

H " Pr Re„ Til M 2 + !- T 

2 00 w 

are plotted in Fig. 9. The present results show reasonable agreement with 
the experimental measurements, although some disagreement is observed 
near the beginning of the ramp. 
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The details of the computed flowfield are illustrated in the contour 
plots of Figs. 10, 11, and 12, Shock details can be seen in Fig. 10 where 
the relatively weak leading edge shock is represented by a single pres- 
sure contour, Pictured in Fig. 11 is the shear layer which emanates from 
the shock intersection as well as the strong density gradient regions in 
the boundary layer and at the induced shock. The Mach contours of Fig. 12 
Illustrate the rapid decrease in the boundary-layer thickness which results 
in the large increase in heat transfer to the plate shown in Fig. 9. It 
should be noted that the vertical coordinate in these contour plots has 

<v 

been unsealed by the factor l. 

In order to Include Influences from downstream and, thus, hopefully 
correct the disagreements shown in Figs. 8 and 9, a global iteration on 
the pressure field was Incorporated into the PNS code. PNS methods have 
proven effective in solving supersonic flows with weak viscous-invlscid 
interactions characterized by small pressure gradients in the marching 
direction. In the present hypersonic test case, however, large pressure 
gradients are generated and the associated upstream influences are there- 
fore expected to be significant. This upstream effect is Ignored in the 
single sweep PNS method and several researchers have turned to global 
iteration techniques to include these effects. 

Recently, Raklch (17) developed an iterative method in which the 
pressure gradient term on the right side of the momentum equation is 
evaluated using a combination of values at two Iteration levels. That 
is, AP in Eq. (lib) would be evaluated with 
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AP 
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(1 - ft» 
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n+2 


n+1 


)■ 


0 , 


0 


T 


J J 

where the subscript i indicates the present iteration level. It is 

necessary here to move the implicit part of the forward difference, p^ + , 

to the left-hand side in Eq. (lib). This is done through a series of 

algebraic manipulations after which the form of the term added to the 

» 

right-hand sides of Eqs. (lib) and (12b) is 


AP = 


0 , - (1 - OJ) 


n+2 n 
-1 v l 


( UTi 
¥ 


T 


0 , 


(17) 


The effect on the left-hand side is to change E* to 

T 

= j [ pu, pu 2 + (2to - l)p, puv , (e fc + p)u j 


The new left-hand sides of Eqs. (lib) and (12b) are then obtained by 
simply replacing W in the A Jacobian with the term (2io - 1) . 

This iterative method is implemented by taking the first sweep with 
the standard PNS code. In subsequent iterations, E* is replaced with E^ 
and AP is evaluated according to Eq. (17), Only the pressure field is 
stored after each iteration. Following Rakich's technique for the outflow 
boundary condition, the pressure gradient is kept constant and equal to 
the gradient calculated at the end of the first sweep. 

The global iteration applied with the implicit MacCormack scheme 
converged very rapidly. Figures 13 and 14 show pressure and heat transfer 
coefficient results of calculations which converged in 6 steps and re- 
quired 70 seconds of CPU time on an NAS AS/6 computer. These figures 
show that the converged results deviate only slightly from the results 
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of the first sweep. Even in the vicinity of the corner where the upstream 
Influence should be strongest, the effect of the ramp seems to be felt 
only one or two grid points upstream. 

For purposes of comparison, the global Iteration procedure was 
applied with the Beam-Warming scheme and calculations of the same test 
case were performed. This method seems to be substantially more sensitive 
to the forward difference on the pressure than the Implicit MacCormack 
scheme. The consequences are a slower convergence rate and a solution 
which better displays the effects of upstream Influences. Figures 15 and 

16 show the pressure and heat transfer coefficients, respectively, after 
25 iterations. Regarding the upstream Influence, the point of Interest 
in Fl^. 15 is the gradual pressure increase in the iterated solution 
beginning several points ahead of the ramp. The associated feature in 
Fig. 16 is the slight dip in the heat transfer coefficient near the cor- 
ner. Both of these characteristics agree well with those of the experi- 
mental results though the calculated dip in the heat transfer is displaced 
slightly toward the ramp. 

The pressure, density, and Mach number contours of Figs. 17, 18, and 
19, respectively, illustrate the details of the computed flowfield after 
25 iterations with the Beam-Warming scheme. These figures, especially 
the pressure contours, reveal some oscillatory behavior also present to 
a lesser degree in the wall coefficient plots of Figs. 15 and 16. However, 
the streamwise oscillations of Figs. 15 and 16 are thought to stem from 
the severity of the initial conditions, whereas the oscillations of Fig. 

17 are probably caused by nonlinear effects associated with the induced 
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shock. These oscillations Illustrate the lack of Inherent artificial 
dissipation In the Beatn-Warming scheme. 

Velocity profiles computed during the first sweep by the Implicit 
MacCormack scheme and the Beam-Warming scheme are shown In Figs, 20 and 
21, respectively. The two figures are similar except for the profile on 
the ramp nearest the corner. For this profile, the Beam-Warming scheme 
results seem to show a slightly greater tendency toward separation than 
the results of the Implicit MacCormack code. An expanded view of the 
velocity profiles in the corner region is given In Figs. 22 and 23. The 
presence of an Inflection point In the Beam-Warming scheme results at 
x/Z = 1.014 Is clear, whereas the results of the implicit MacCormack 
scheme remain nearly linear throughout this region. Thus, the code using 
the Beam-Warming scheme, even in the first sweep, seems to exhibit a 
greater sensitivity to the sharp corner than the implicit MacCormack 
scheme code. 

An iterated Beam-Warming scheme calculation produced the results 
shown In Fig. 24. Comparison with Fig. 23 shows that the effects of 
iteration are, again, to propagate the effect of the ramp upstream and, 
in this case, to slow the flow near the wall. The results of the iterated 
implicit MacCormack scheme calculation were indistinguishable from those 
of Fig. 22 and, therefore, have been omitted. 

One source of error in the implicit MacCormack code may, as In the 
flat plate boundary-layer test case, be the solid wall boundary treatment. 
For that test case, the reflective condition, Eq. (16), proved to be the 
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more stable and accurate procedure, However, in the present test case, 

application of Eq. (16) at the wall seems to trigger instabilities near 

the beginning of the calculations. For this reason, the quantity 
f I | — ] n+1 

|B|6Uj^ was set to zero in these computations. Thus, the boundaries 
are treated only explicitly in the implicit MacCormack code whereas the 
implicit conditions used in the Beam-Warming code caused no difficulty for 
this case. 

Another possible source of error in the code employing the implicit 
MacCormack scheme is the effect of subjecting the governing equations in 
strong-conservation-law form to a finite-difference calculation without 
properly satisfying the geometric conservation law. This potential prob- 
lem has been discussed previously in relation to Test Case I where small 
amounts of grid stretching caused very large errors in that high Courant 
number case. In the present test case, however, a Courant number of 
approximately 20 was the largest encountered and only moderate clustering 
of the grid was employed. In addition, "differencing the freestream" 
seemed to indicate that these effects should not be significant. Never- 
theless, it is possible that the nonlinearity of the solution itself 
would accentuate these nonphysical effects and cause observable errors in 


the solution. 
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V. CONCLUDING REMARKS 

The Implicit MacCormack scheme has been applied to the two-dimen- 
sional, parabolized Navier-Stokes equations for the computation of 
steady supersonic flowfields. In order to test this method, two flow- 
fields were computed. These Included a laminar flat plate boundary-layer 
case and the hypersonic laminar flow over a 15° compression corner. Pres' 
ent results compare very well with previously published computational 
results and experimental data. Comparisons were also made with results 
of the conventional Beam-Warming scheme in terms of accuracy, stability, 
computer time, computer storage, and ease of implementation. 

Also presented in this report are results of a global iteration 
technique which allows upstream influences to be included in PNS calcula- 
tions. Very good results were obtained when this procedure was imple- 
mented with the Beam-Warming scheme, but the response of the implicit 
MacCormack scheme to the iteration was somewhat disappointing. In order 
to improve the accuracy of the implicit MncCormack scheme in general, 
further investigation is suggested into the development of a reliable 
implicit boundary condition procedure. In addition, to reduce the strong 
grid dependence observed here, use of either the chain-rule-conservation- 
law form of the equations (see Refs. 16 and 20), or a finite-volume 
approach (see Ref. 13) is strongly recommended. 
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where $ = -r(y-l) (u + v ). Defining V = 11 u + T) v, 9F./91J can he written 



the matrix of left eigenvectors of 3F./3TJ is 
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